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1. OBJECTIVE 

The objective of this research is to control the dynamics of multizone 
programmable crystal growth furnaces. Due to the inevitable heat exchange among 
different heating zones and the transient nature of the process, the dynamics of 
multizone furnaces is time varying, distributed and therefore complex in nature. 

Electrical power to heating zones and the translational speed of the ampoule are 
employed as inputs to control the dynamics. Structural properties of the crystal is the 
ultimate aim of this adaptive control system. These properties can be monitored in 
different ways. Following an order of complexity, these may include : 

a. - on line measurement of the material optical properties such as the 
refractive index of crystal, 

b. - on line X-ray imaging of the interface topology, 

c. - on line optical quantification of the interface profile such as the 
determination of concavity or convexity of the interface shape, and finally, 

d. - on line temperature measurement at points closest to the material such as 
measurements of the ampoule's outside and inside surface temperatures. 

The research performed makes use of the temperature and optical 
measurements, specified in (c) and (d) above as the outputs of furnace dynamics. 
However, if the instrumentation is available, the proposed control methodology can be 
extended to the measurements listed in (a) and (b). 

2. ACCOMPLISHED RESULTS 

Following results are accomplished in this study. 

A - Design, construction and validation of a computerized 8 zone heater drive 
and monitor system. The system has the following technical specifications. 

Power drive capacity : 2400 W 

Temperature reading resolution : 16 bits for 0 to 

1500 C temperature range. This corresponds to 0.02 C 

temperature reading resolution. 

Heater drive resolution : 16 bits. 

Temperature accuracy at the highest setting : 1 C. 

Optical resolution of the interface detection : 12 bits 

B - Experimental identification of multivariable furnace dynamics. 
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Pseudo random binary heating power is introduced to the furnace. The 
multivariable transfer function of the furnace is determined between the input and the 
resulting zone temperatures. 

C - Based on the identified model structure, design and verification of a 
multivariable controller is accomplished. 

3. INTRODUCTION 

One common way to grow crystal is to subject the molten material to a 
temperature gradient. This requires keeping a predetermined axial and radial 
temperature profile around the material despite the effects of unknown process 
disturbances. To maintain continuous solidification, either the material or the 
temperature gradient is moved axially. At any given point, the temperature is a function 
of axial position h, radial position r, circumferential position c, and the time t. This 
distributed nature of the temperature can be described by a multivariable function 

T = T(h,r,c,t) 

The time dependency of the temperature is due to the following factors: 

i- the translational motion of the ampoule with respect to furnace, 

ii- the inevitable heat exchange between different heating zones. 

The temperature control problem can now be described as follows. Determine 
the electrical heating power to individual heating zones and the translational motion 
velocity of the ampoule such that at each point within the ampoule the temperature, 
T(h,r,c,t), is kept at the desired level, T d (h,r,c,t) ,i.e, 

T(h,r,c,t) = T d (h,r,c,t) 

Previous research has shown that the shape of the solid liquid interface is 
strongly correlated to the imposed temperature gradient and the translational motion 
velocity of the ampoule. Furthermore, the macroscopic shape of the interface affects 
the crystalline perfection of the material. It is known that [Chang and Wilcox 1974] 
minimization of thermal stresses can be achieved by maintaining a planar interface 
between the solid and the liquid phases of the material. Such planar interface is 
desirable because crystalline imperfections are also minimized. [Taghavi and Duval 
1989] have studied the effects of furnace temperature profile and translational velocity 
on the interface shape, and they have also considered the inverse problem of the 
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desired temperature profile that would result in a planar interface. These studies have 
given insight into the problem of achieving planar interface during crystal growth . 
Most importantly, the ability to design a furnace where the temperature profile can be 
tailored during the growth of the crystal is found to be highly desirable for obtaining a 
planar interface. This has led into the concept of transparent multi-zone Bridgman 
furnace. 

Experimental observations confirm that the same factors affect the shape of 
interface, namely : 

a. the axial temperature profile alongside the ampoule, 

b. the radial temperature profile surrounding the interface, and 

c. the translational velocity of the ampoule. 

[Wang et al 1984] experimented with a single zone vertical Bridgman type 
furnace using Gallium doped Germanium as the material. They considered the energy 
flow, the temperature distribution and the translational speed as their main parameters. 
They verified that the macroscopic growth rate was always transient and the interface 
remained concave with respect to solid at the translational speed of 4 micro meter/sec. 
They further concluded that the main parameters effecting the crystal growth 
process are : 

i) the establishment of a well defined and symmetrical temperature profile 
around the ampoule, and 

ii) the provision for the growth interface quantification and control. 

The Mellen EDG furnace [Parsey and Thiel 1985] can emulate the heat flow 
characteristics of the horizontal Bridgman-Stockbarger furnace without requiring the 
motion of the ampoule. The use of the Mellen furnace in horizontal configuration for a 
large diameter Gallium Arsenide crystal has been investigated by [Bourret et al 1987], 
They have found that the crystalline perfection is improved by increasing the axial 
temperature gradient over the solid and decreasing the radial temperature gradient over 
the melt. 

Quantification of solid-liquid interface by image processing has been performed 
by [Batur et al 1990a, 1990b]. This powerful technique can locate the position of 
interface and determine its concavity or convexity for transparent furnace and 
materials. This topological information about the interface can serve the following 
purpose. 

- It provides continuous on line information for the position and the shape of the 
interface. Therefore the effects of different temperature gradients and the translational 
speeds on the interface can be quantified in real time. 



6 


- It can be used as the feedback signal in conjunction with the classical 
temperature feedback from each zone. 

In conclusion, material structural properties, partially observed by the shape of 
solid-liquid interface, are complex, time dependent phenomena of the temperature 
profiles and the ampoule's translational velocity. A multivariable experimental furnace 
model explaining the interactions between different heating zones and the effect of the 
translational velocity is the first fundamental step to understand the process. The 
purpose of the proposed investigation is to find this model in real time and based on the 
model, implement an adaptive multivariable control algorithm for temperatures and 
interface profile. 

4. CRYSTAL GROWTH CONTROL 

Control of crystal growth should ideally be based on the desired structural 
properties of a given material. For example, if the acousto-optical properties of the 
material are critical then the control law should be designed around these parameters. 
On line measurement of the solid's refractive index, for example, provides such a 
feedback signal for the control algorithm. However if the desired structural properties 
can not be directly measured because of cost or furnace design, then two alternative 
options exist. If the relation between the desired state and a physically realizable 
measurement is known then this relation can be used to find the desired measurement. 
If this relation is unknown, the next best measurement, closest to the desired 
observation should be provided for the controller. 

The shape of the interface is one measurement which is closely related to the 
structural properties of the material. X-ray imaging of the ampoule, for example, may 
provide information to determine the solid-liquid interface. The differences in X-ray 
intensities corresponding to solid and melt can be exploited to determine the topology 
of the interface. X-ray imaging is an elaborate and equipment intensive process. For 
transparent furnace and material, if there is significant difference in gray levels or 
colors between the solid and the melt, optical imaging can also locate the position and 
the shape of the interface [ Batur et al 1990a, 1990b], 

Perhaps the most practical furnace control methodology is to measure and 
control the temperatures which define the necessary thermal gradient for the growth of 
the crystal. Although this control technique is not based on the measurement of the 
desired structural properties of the material, nevertheless it is preferred for the 
following reasons. 
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- As outlined in the introduction, analytical studies have shown that the shape of 
the interface can be related to the axial and radial temperature distribution surrounding 
the interface. 

- Temperature measurement technology is well established and reliable in 
comparison with, for example, optical measurements related to the crystal properties. 

In the next section we will describe temperature based control methodologies for 
the crystal growth furnace. 

5. CONTROL ALGORITHM 
Problem Formulation 

The system to be controlled is assumed to have the following model structure 
A(z _1 )y(t) = B(z _1 )u(t) + C(z'l)n(t) (1) 

where 

A(z'l) = I+A 1 z‘ 1 +A 2 z- 2 + ....+A na z- na , Aj e RP X P 

B(z-!) = B 0 z- k +B 1 z' k - 1 + +B nb z- k * nb , B[ eRP x P 

C(z _1 ) = I + Cjz' 1 +C 2 Z' 2 + + C nc z’ nc , Cj e RP X P 

where u(t), y(t) and n(t) are p x 1 vectors representing input, output and the disturbance 
of the system, k is the process dead time. The polynomial A(z'l), B(z"l), C(z 1) may 
be time varying. 

The resulting control law takes the following form. 

R(z-l)u(t) = T^Mt) - S(z*l)y(t) (2) 

where R(z _1 ), T(z _1 ), S(z _1 ) matrix polynomials are designed in such a way to achieve 
desired characteristics as discussed later and they are assumed to take the following 
form 

RCz- 1 ) = I+R 1 z- 1 +R 2 Z- 2 + .... + R nr z- nr , Rj e RP X P 
T(z _1 ) = To+T 1 z- 1 +T 2 z-2 + ....+T nt z- nt , Tj e RP X P 
SCz' 1 ) = So+S 1 z- 1 +S 2 Z- 2 + ....+S ns z- ns , Sj e RP X P 
r(t) is the desired reference level of the controller. 

Given the system as in (1), it is desired to find a self-tuning control law as in 
equation (2), so as to make the input sequence u(t) and the tracking error e(t) defined as 
[r(t) - y(t)] as stable sequences, independent of C(z~l)n(k) and r(k). It is also assumed 
that A, B, C may be slowly time varying. 

Controller Design Problem 

The closed loop equation is obtained using equation (1) and (2) as 
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Ay = BR _1 (Tr - Sy) + Cn 

(A + BR^S) y = BR^Tr + Cn (3) 

The dynamics of the closed loop system is dictated by (A + BR _1 S) matrix 
polynomial. This can be made equal to any desired polynomial by suitably choosing R 
and S. Unfortunately, in trying to obtain R and S using the above method, we end up 
with a set of non-linear equations and there are no standard techniques to solve this 
type of equation. Hence, the following transformation is performed 

R-iS = SR' 1 (4) 

[Wolovich 1974] has shown that equation (4) has an unique solution if R and S are 
given as 

R = I+R 1 z' 1 + ....+R nr z' nr , Rj e RP X P 
S = I+S 1 z' 1 + .... + S nr z' nr , Si e RP X P 
From equation (4), it follows 



Hence 

det(R) = det(R) (6) 

because R and R have the identity matrix as its first coefficient. From (4) and (6), it 
follows 

det(S) = det(S) (7) 

Using transformation (4) in (3), we obtain 
(A + BSR'^y = BR _1 Tr + Cn 

( AR + BS) R'iy = BR' J Tr + Cn (8) 

In the pole placement procedure, R and R are found to satisfy the following 
design equations 

AR -1- BS = P (9) 

where P is the desired polynomial matrix. As observed in (8), the dynamics of the 
closed loop system is also described by R' 1 . In the later section, it is shown that this 
does not affect the stability of the dosed loop system, even though R itself may be 
unstable. 

The main idea behind this algorithm is to solve for R and S matrices in equation 
(9) and then use these to find R and S using (4). Before getting into details of the 
solution to equations (4) and (9), let us consider the steady state performance of the 
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regulator. In order to have no steady state error, R matrix is forced to have an 
integrator in each loop as given by 

Rfz* 1 ) = (l-z'l)R'(z' 1 ) (10) 

As mentioned before, we have no direct method to select R matrix and its value 
is determined through R. Fortunately equation (6) comes to the rescue, which 
essentially states that the roots of characteristic equations of R and R matrices are the 
same. Hence R also has a forced integrator and is given as 

R(z~l) = (l-z'l)R'(z'l) (11) 

with 

R' = I+Rjz- 1 + ....+Rn r -iz- nr - 1 , Rj e RP X P 
Now the steady state conditions for step inputs in r(t) is obtained from equation (8) by 
considering the limit z - > 1. From equation (11), it follows that R tends to zero when 


z - > 1 . Hence the steady state conditions are given as 

[ A(1)R(1)+B(1)S(1)] R'*(l)y(t) = B(l)R l(l)T(l)r(t) (12) 

B(l)S(l)R-l(l)y(t) = B(l)R-l(l)T(l)f(t) (13) 

Applying transformation in equation (4) to (13), we get 

B(l)R-l(l)S(l)y(t) = B(l)R-l(l)T(l)f(t) (14) 

Hence 

y(t) = S-l(l)T(l)f(t) (15) 

Now the steady state tracking error is given as 

e(t) = f(t) - y(t) = [I - S-l(l)T(l)]f(t) (16) 

e(t) = S'l(l)[ S(l) - T(l)]f(t) (17) 

Again T matrix is selected such that 

T(l) = S(l) (18) 

As a result of the above constraint on T 

e(t) = 0 (19) 

Instead of just satisfying equation (18), we pick T as 

T(z-l) = S(z-l) (20) 


Although this ensures zero steady state error for step inputs, the excess freedom 
in choosing T matrix is not properly utilized. [Mikles 1990] exploited this to follow 
ramp and other higher order inputs. At this point, it is just noted that this technique 
can be modified for higher order reference tracking problems. 


Solution to Design Equations 



10 


R and S are found by solving (9) with R having the form as in (11). 
Substituting (1 1) in (9), we get 

AR'O-z* 1 ) + BS = P (21) 

This can be solved by using Toeplitz form to represent the above relation. 
Introducing 


r Oi 


r Oi " 

0 V 


Ok 

I 


BO 

Ai-I 


Bl 

J A 2 -Ai 


b 2 K 

Ana' Ana- 1 


Bnb-1 

'Ana 


®nb 

Ona+v+2 


Onb+k+v 

" 0 nn J 


“ o nn > 


where np is the order of P matrix and each Oi is a pxp null matrix. 


( 22 ) 


Then the 


system denoted in (21) is given by 


[ A0....A nr _iBQ....B ns ] 

r 2 


r p > a 
?2 

Al' 1 \ 
A 2 -Ai 


R nr 


*• 

Ana* Ana- 1 


so 


•' 

'Ana 


Si 


•* 

Ona+2 




R np-1 

... 


V^ns/ 

V R np / 

V o n p J 


This is a linear system of equation and solutions to such equations are well 
known. In particular, LU decomposition is used to solve for the unknown matrices. 

In a similar fashion equation (4) is solved with R having a form as (10) and R 
satisfying (1 1). Substituting relations (10) and (1 1) in (4), the following result is 
obtained. 

(l-z- 1 )- 1 R*- 1 S = SR’-^I-z' 1 )* 1 


( 24 ) 
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Cancelling the common terms in (24), we obtain 

R'^S = SR'' 1 (25) 

It can be seen that (25) is similar to (4). Rearranging and then taking transpose, 
equation (25) takes the following form 


r'T s T . sTr’T = o (26) 

Here R*, S are known quantities and R , S are the unknown quantities for which 
the equation has to be solved. This equation is similar to (21) and hence adapting a 
similar procedure, we define the following matrices. 


W 


r Oj ^ 


' o, ' 

0 V 



I 

- .T 


?¥ 

<j R 1 

► Sv = < 

Si ► 

-V T 


-T 

R nr 


Sns 

°nr+v+l 


Ons+v+1 

Onr+ns 


^ Onr+ns 


(27) 


with each Oj being a pxp null matrix. Using the above relation (26) can be represented 


in a Toeplitz form as 

T T T 

[ R 0 •••• R„s'Si 


T 

-Sir! 


S 1 


R 1 


T 

= [S }] 


(28) 


By solving this equation, we obtain the control law polynomials R and S. 

Self Tuning Regulator 

The off-line version described in the previous section assumed the knowledge of 
A, B, C polynomials. These are not known off hand but can be estimated. There are 
many algorithms available to estimate such systems depending on the type of the C 
polynomial. In our application, i.e. identification of dynamics in the furnace, it is 



assumed that C=I and hence recursive least squares algorithm is used to identify the 
system parameters. 

Recursive Least Squares Algorithm 

The system is assumed to have a model structure as in equation (1) 


A(z _1 ) y(t) = B (z _1 )u(t) + C(z-!)n(t) (29) 

with all the symbols same as before. Converting (29) into the vector difference 
equation notation as 

y(t) = 'A iy(t-l )-.... -A na y(t-na) 


+ B ! u(t- 1 ) + . . . + B nb u(t-nb) + n(t) (30) 


Defining the following matrices 


Y v = 

" y T (-v) ~ 
y T (l-v) 

u v = 

uY(-v-k) " 
yY(l-v-k) 

V = 

~n T (l)~ 

n T (2) 


-y T (N-v)- 


-y T (N-v-k)- 


-nT(N)_ 


Where N is the total number of observations. Using notations in (31), y(t) for N 
observation is given as 


[-Y 1 -Y 2 ...-Y na U 0 U 1 ...U nb ] 



Y 0 = 



+ V 


(32) 



Defining 

f= C-Yj -Y 2 - -Y na U 0 Ui ... U nb ] 
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Complete Self-tuning Algorithm 

1. Assume initial qQ, P n =al, where a is a large number. We are not sure about 
the initial value of q and so by keeping P large, we have a high initial gain in 
the updation equation which result in the estimated parameters having a very 
fast initial convergence rate. Actually large a implies very little confidence on 
the initial guess of qg. 

2. Past value of y and u are used to update q and P using equations (38) and (37). 

3. Using q, we find A(z'l) and B(z‘l) from (33). 

4. With this A(z _1 ) and B(z _1 ), R' and S are found with the help of (28). 

5. Now r' and S are found using (28). 



6. R* and S are the necessary matrices needed to determine the control input u(t) 

using (2). This input is applied to the system. 

7. Repeat from step 2. 

The self-tuning property of this algorithm has been proved by [Prager and 
Wellstead 1980]. They assumed the convergence of the system parameters to show this 
feature of the algorithm. The overall convergence itself is yet a subject of research 
even in SISO models. But many authors have shown this algorithm working on 
simulation and there are few practical application which have made use of this 
algorithm. 

Stability of overall closed loop system 

In the pole placement problem formulation, we have stated that the designed 
controller must generate a control input u(t) and the tracking error e(t) as stable 
sequences. In the previous section, we derived a self tuning controller based on the 
pole placement design. In this section, we will show that u(t) and e(t) are stable 
sequences. 

The system , controller, transformations and pole placement design equation are 
given below for convenience with all the notations same as before. 


System 

Ay = Bu + Cn 

(39) 

Controller 

Ru = Tr - Sy 

(40) 

Transformation 

R-iS = SR' 1 

(41) 

Design Equation 

AR + BS = P 

(42) 

Introducing additional transformation equations 


A^B = BA* 1 


(43) 

RA + SB = P 


(44) 



A has 


similar structure as R with the first term being the identity matrix and has an order n a . 
Hence, in a similar fashion as in (5) and (6), the following relations can be derived. 


det(A) = det(A) (45) 


det(B) = det(B) 


(46) 
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Now equation (42) can be written as 
P = AR + BS 

= bib-Ur + S) ( 47 ) 

Using the fact in (43), (47) can be rewritten as 
P = B(AB _1 R + S) 

= B(AB' 1 + SR _1 )R 

Again using the transformation in (41), we get 
P = B(AB"1 + R _1 S)R 
= BR‘ 1 (RA + SB)B*^R 
From (44), we get 


P = BR'lPB'lR ( 49 ) 

Hence 

det(P) = det(B)det(R' 1 )det(P)det(B' 1 )det(R) (50) 

From relations in (5), (6), (45), (46) and using the fact det(R ^) = 

det(P) = det(P) ( 51 ) 

Equation (39) can be written as 

y = A'^Bu + A'^Cn (52) 

Substituting (52) in (40) and multiplying the resulting expression by R _1 , we get 

u = -R-iSA-iBu + R _1 Tr - R^SA^Cn (53) 

Further Simplifying (53), we get 

u = (Im+R-lSA-iBr^Tr - (Im+R-iSA^B^R^SA^Cn (54) 

Now (I m + R-1 sA‘1b)" 1R"1 can be simplified as 

(I m + R-lSA-lBrlR' 1 = [R + SA^B]- 1 (55) 

Taking advantage of transformation defined in (43), we get 
(I m + R'iSA^B)-^- 1 = [R + SBA' 1 ]’ 1 
= ([RA + SB] A’ 1 )' 1 

= A (RA + SB)' 1 ( 56 ) 


In a similar fashion (I m + R^SA^B^R^SA' 1 can be written as 
(I m + R^SA^B^R^SA- 1 = [AS _1 R + B]' 1 
= [ARS' 1 + B]' 1 
= ([AR + BS] S' 1 )’ 1 
= S (AR + BS)' 1 

Substituting (56) and (57) in (54), u becomes 


( 57 ) 
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u = A (RA + SB)'^Tr - S (AR + BS^Cn (58) 

From the relation in (42) and (52), u can be written as 

u =AP- 1 Tr-SP- 1 Cn (59) 

Since P and P are designed and are selected as stable matrices; the resulting 
input sequence is stable. 

Now the tracking error e is given by 
e = r-y 

Using the relation in (44), we get 

e = r - A _1 Bu - A^Cn (60) 

Substituting the expression (58) for u in the above relation, we get 
e = r - A-!b( A (RA + SB^Tr - S (AR + BS^Cn) - A^Cn 
= (I r - A-lBA (RA + SB)'^T)r 

+ ( A _1 B S (AR + BS)" 1 - A' 1 ) Cn (61) 

The first term can further be simplified as follows 
(I r - A* 1 BA (RA+SB)- 1 T) 

= (I r - BA-U (RA+SB)-!t) 

= I r -BP- 1 T (62) 

Similarly, the second term in equation (61) is simplified as 
(A^B S (AR + BS)- 1 - A' 1 ) 

= ( A'iB S - A- ] (AR + BS)) (AR + BS)' 1 
= -R(AR + BS)’ 1 

= -RP’ 1 (63) 

Hence the tracking error e is given by 

= (Ir-BP^TJr-RP^Cn (64) 

Again since P and P are designed as stable matrices, the resulting tracking error 
sequence is stable. Hence, under closed loop, both the input and the error sequences 
are stable. 


Stability of Controller 

In the previous section, we proved that the overall closed loop system is stable. 
There are some concerns about the controller itself becoming unstable. The input is 
given by 


u = R -1 Tr - R -1 Sy 


(65) 
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R already has an integrator which corresponds to a pole at | z | = 1 . The other 
poles of R may become unstable i.e. \z\ > 1. In other words, the poles of R can be 
anything as there is no provision in the design that makes R necessarily stable. This is 
of big concern in spite of the proof on u being stable sequence irrespective of what R 
is. The proof assumed zero modelling errors, unlimited inputs. In practice, we always 
have limited inputs and especially during self tuning mode, we have modelling errors. 

When saturation of the input occurs, the system is essentially operating under 
open loop and the resulting system may operate in an undesired manner. Even when 
all the poles of R are stable, the pole at z = 1 can cause high frequency oscillations due 
to a phenomenon known as integration windup. Although many practical solutions for 
windup are known, the presence of unstable poles in R can render these solutions 
ineffective. Some of the techniques to overcome windup are given in the next section. 
When the controller input gets saturated, the actual input applied to the system does not 
belong to the controller. If after sometime, the input do return to the unsaturated 
region, the past inputs and outputs in equation (40) do not belong to the controller itself 
as a result of which the output of the system can take undesirable tours. 

At this point, we conclude that if all the poles of R are at least marginally 
stable, then we can take care of limited inputs as given in the next section. 

Some Practical Solutions to Integration Windup 

Integration windup occurs, when the control input demanded by the controller 
saturates the actuator. The easiest solution is to remove the integrator from the control 
loop whenever saturation in the actuator occurs. This can be done by directly solving 
(19) without additional constraints given in equation (10). In this situation, we have 
two controllers with only one controller operating at any instance. 

Another scheme which is more attractive is to include another loop in the 
control system as given in Figure 2. This loop gets activated only when the input 
demanded by controller u a . The individual loop itself can be seen in Figure 3. The 
transfer function between u c and u a is given from Figure 3. 

u s = u a - u a (66) 

and u c = R-^FCua - u c ) + z) (67) 

Hence, the transfer function between u a and u c is given by 

u c = (R + F)-lFu a (68) 

The main aim of this loop is to make u c -> u a and this is done by suitably 
designing F so as to make (R + F) asymptotically stable. This section gives us some 
techniques to overcome integration windup. The effect of the new loop on the overall 
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stability of the closed loop system is not shown here as the purpose of this section is to 
list certain techniques to get around integration windup. Again, if R has unstable 
poles, it is quite possible that these techniques may prove ineffective. 


Example 

In this section, examples are given to illustrate the pole placement design 
procedure in diagonal and non-diagonal systems. 

(a) Diagonal Case 

The diagonal system is given below 


1-1.5z"1+0.7z'2 

0 


l-1.5z'l+0.7z' 2 J y ^ . 


r 1 +z' 2 
0 ; 


0 

:‘l +z‘ 2 _ 


u(t) 


(69) 

We want to place all the poles of the system at z= +0.5. In order to achieve 
this, we have to find R’ and S satisfying the design equation given in (21). R . and S 
are assumed to be second order polynomials. Hence P, the desired polynomial must be 
of 5th order as given below 

J (l-0.5rl)5 0 I (70) 

v L 0 (l-0.5z _1 )^J 

Using the procedure described in the earlier section, R and S are given as 
_ ri-0.1419z* 1 +0.04464z- 2 0 1 

“L 0 1-0. 1419z -1 +0.04464 z‘2j 

0.142-0.241z- 1 +0.115z- 2 0 1 

0 0. 142-0. 241z _1 +0.1 15z' 2 J 

(71) 

(b) Non-diagonal Case 

The system is given by 


and 


R 

S 


= [ ( 


' l-0.5z _1 

1+0.5Z- 1 ' 

y(t) = 

z-i 

0.2z _1 " 

. 1-0. 5Z' 1 

l-0.5z*l . 

i 

O 

LP 

N 

i 

►— * 

Z“^ . 


u(t) 


(72) 


Again, we want to place all the poles of the system at z= +0.5. R and S are 
assumed to be first order polynomials. Hence P, the desired polynomial must be of 

3rd order as given below 
(1-0.5Z- 1 ) 3 


P = 


0 


0 (1-0.5Z- 1 ) 3 J 


( 73 ) 


R* and S are found using equation (21) and 
_ [ 1-0.729Z’ 1 1+0.729Z' 1 ' 

R “ U-0.729Z- 1 1-0.729Z’ 1 . 
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and S = 


(74) 


-0. 574+0. 575z'l -0.093+0.093Z- 1 
L 0. 344-0. 344Z' 1 -0.452+0.453Z- 1 J 
Simulation of this non-diagonal system for tracking references for two different 
pole location 0.5 and 0.95 are given in Figures 4 and 5. While placing poles at 0.5, 
the resulting controller is stable and with 0.95, the controller is unstable. This 
controllers are made to track references, some of which cannot be achieved with the 
limited inputs. It can be seen from the figures that the stable controller has much more 
desirable characteristic than the unstable controller. There are lot of oscillations when 
saturation in the actuators occurs in the case of unstable controller. For this reason 
stable controllers are preferred. 

In our application domain, there are 8 inputs and 8 outputs. It is also observed 
that each term is at most a 2nd order polynomial. Hence, the furnace is modelled as a 
8x8 MIMO model as given below 

y(t) = -A 1 y(t-l)-A 2 y(t-2)+B 1 u(t-l)+B 2 u(t-2)+n(t) (75) 

where y(t) and u(t) e R^ x l and A], A 2 , Bj, B 2 e R^ x ^. Here the elements of matrix 
R and S are assumed to be at most 2nd order polynomials, ie 
R' = I + Rjz- 1 + R 2 z- 2 

and S = Sq + Sjz"^ + S 2 z"2 (76) 

The desired polynomial matrix P is a diagonal matrix with all diagonal terms 
being 5th order polynomials. This is used to control a simulated furnace. Preliminary 
actual input-output data of the furnace under self tuning pole placement control is given 
in Figures 6 and 7. It can be observed that initially there is lot of fluctuation in the 
input and the input hitting the limits. But after some time when the estimated 
parameters converge the oscillations are reduced and the desired temperature level is 
achieved. It is also observed that there are some oscillations about the reference 
temperature. It is believed that these are caused due to poor estimates resulting from 
the poor initial guess on the system parameter. 


6. APPLICATION TO CRYSTAL GROWTH FURNACE 


Application of the control technique developed in this research is presented in 
Appendix 1. 

7. CONCLUSION 

In this reserach, a self-tuning pole placement control algorithm is developed for 
the temperature control problem in crystal growth furnaces. The algorithm is applied 
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successfully to control the axial temperature around the ampoule. Its suitability in 
controlling the temperature of a 8-zone crystal growth furnace is discussed. The 
stability of the overall control system is derived. Some of the robustness issues 
concerned with this algorithm are analyzed. 
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This paper presents an intelligent adaptive control system for the control of a solid-liquid interface of a crystal while it is 
growing via directional solidification inside a multizone transparent furnace. The task of the process controller is to establish a 
user-specified axial temperature profile and to maintain a desirable interface shape. Both single-input-single-output and 
multi-input-multi-output adaptive pole placement algorithms have been used to control the temperature. We also describe an 
intelligent measurement system to assess the shape of the crystal while it is growing inside a multizone transparent furnace. A color 
\ideo imaging system observes the crystal in real time, and determines the position and the shape of the interface. This information 
is used to evaluate the crystal growth rate, and to analyze the effects of translational velocity and temperature profiles on the shape 
of the interface. Creation of this knowledge base is the first step to incorporate image processing into furnace control. 


1. Introduction 

One way to produce a crystal is to subject the 
pure molten compound, encapsulated inside an 
ampoule, to a temperature profile which includes 
its solidification temperature. This requires main- 
taining a predetermined axial and radial temper- 
ature profile around the material, despite the 
effects of unknown process disturbances. To 
maintain continuous solidification, either the am- 
poule or the temperature gradient is moved along 
the longitudinal axis of the furnace. 

Singh et al. [1] have shown that the shape of 
the solid-liquid interface is determined by the 
magnitude of the imposed temperature gradient 
and the translational velocity of the ampoule. 
Minimization of thermal stresses can be achieved 


by maintaining a planar interface between the 
solid and liquid phases of the material (Chang 
and Wilcox [2]). Also, a planar interface can be 
desirable because of minimized crystalline imper- 
fections. Recently, Taghavi and Duval [3] have 
studied the effects of furnace temperature profile 
and translational velocity on the interface shape. 
They also considered the inverse problem of the 
desired temperature profile that would result in a 
planar interface. 

Wang et al. [4] experimented with a vertical 
Bridgman-type furnace where they studied the 
growth of gallium-doped germanium. They veri- 
fied that the macroscopic growth rate was always 
transient and the interface remained concave with 
respect to solid at the translational speed of 4 
mm/s. 
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The Mellen EDG furnace (Parsey and Thiel 
[4]) can emulate the heat flow characteristics of 
the horizontal Bridgman-Stockbarger furnace 
without requiring the motion of the ampoule. 
The use of the Mellen furnace in horizontal con- 
figuration for gallium arsenide crystal has been 
investigated by Bourret et al. [6], They have found 
that the crystalline perfection is improved by in- 
creasing the axial temperature gradient over the 
solid and decreasing the radial temperature gra- 
dient over the melt. 

Since the shape of the interface is a significant 
characteristic of this process, automatic detection 
and quantification of this variable can play an 
important role in the crystal growth control. Some 
of the advantages include: 

U) The quantification process provides on-line 
information about the position and the shape of 
the interface. Therefore, the effects of different 
temperature profiles and translational speeds on 
the interface shape can be determined in real 
time. 

(2) The crystal growth rate can be calculated by 
detecting the position of the interface. 

(3) Most importantly, the interface shape infor- 
mation may be used as a feedback signal in the 
interface controller. 

In summary, material structural properties are 
complex, time dependent phenomena of the tem- 
perature profiles and the ampoule’s translational 
velocity. The purpose of this study is to find a 
real time model of the furnace and based on this 
model, to implement adaptive control algorithms 
for. the interface control. 


2. Identification of furnace dynamics 

2.1. Identification of single-input -single-output 
(SI SO) models 

The following time invariant, thermal model is 
assumed: 

y,(r)=B,(z- 1 ) «,-(/ - 1) + «,(')> 0) 

where u, )\ and n are the energy inputs, temper- 
ature, and disturbance terms respectively. 


and B,(z~ ] ) are the polynomials of (r*" 1 ) with 

coefficients [a ] a 2 an ^ [^o ^1 

They define the* dynamics of each zone. The 
subscript i refers to a particular heating zone. 
The z' 1 is the delay operator. The disturbance 
a .( r ) is assumed to be a zero mean stationary' 
stochastic process. 

In the furnace studied, motion of the ampoule 
with respect to furnace changes the dynamics of 
each zone. This is due to the changes in the 
thermal capacitance in each zone. Because the 
ampoule's motion is slow. 2-50 cm/day, a self- 
tuning temperature control algorithm is em- 
ployed. As shown below, the process parameters 
for each zone are estimated by an on-line sequen- 
tial least squares algorithm (Goodwin [7]): 

0,(0 -fli('-i) 

+ 1 + - 1) p,(l ~ 1) d>,( t - !) 

x[y,(0-e, T (r-i) <^,('- 1 )]. ( 2 ) 

p,( 0 

-P,0 -1) 

P,(t - 1) Sfit - 1) - 1) Pfj - i) 

l + *7(r-l) <*>,(' -n 

(3) 

e, = [flj a 2 ••• a na b 0 b\ ■■■ ( 4 ) 

(f),.(/-i) = [-.v(J-i) -y(> - 2 ) • 

-y(r - na) u(t - 1) u(t - 2) 

X • • ■ u(t -nb)]] , (5) 

where 0, is the process parameter estimates and 
( p t (t — 1) is the data vector for zone /. The initial 
choices for the algorithm are 0(0) - 0 and P( 0) = 
al, where a is a large number indicating little 
confidence on the initial estimates. 

2.2 Identification of multi-input -multi-output 
(MIMO) models 

The single-input-single-output (S1SO) models 
do not explicitly account for the temperatures 
and input powers of the neighboring heating 
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zones. Changes in power and temperature of 
adjacent zones act as inputs to the controlled 
zone. Any one zone may be affected by the power 
and temperature of other zones, which may act as 
a disturbance. This phenomenon is observed 
through a non-proportional power distribution 
along the controlled furnace. The power distribu- 
tion affects the amount of radiant energy trans- 
mitted into the crystal from a particular zone. For 
these reasons, a multivariable model has also 
been determined for the system. 

A multivariable model takes into account the 
interaction between the zones. By knowing the 
interactions, the controller can be designed to 
further minimize the temperature fluctuations 
within the furnace. The input and output dynam- 
ics of the system are represented by the following 
multivariable model: 

A(z~') y(r )-B(z-‘) u(r-l)+#.(r), (6) 

where A(z~ ] ) and B(z~ l ) are polynomial matri- 
ces as, 

na 

A(z~')=I+ (7) 

/«! 

nb 

B( 2 ->)= Es,z-', (8) 

i-0 

where A, and B i in (7) and (8) are all 8x8 
matrices. For an eight-zone furnace, input-out- 
put and noise vectors are 


y (0 = [>’.( 0 y 2 (0 •• 

• Va ( 0 ] T . 

( 9 ) 

“(0 = [ w i (0 M :(0 ‘ ‘ 

' “s(0lk 

( 10 ) 

«(0 = [« i (0 «:(0 ■ 

’• «8(0] T - 

( 11 ) 


Note that this multivariable model is drasti- 
cally different from that of the single input single 
output. Zone interactions are explicitly modeled 
in (6) as off diagonal terms of the A x and B t 
matrices. The S1SO model, on the other hand, 
tries to accommodate these interactions by addi- 
tional noise term in the model 

The multivariable model is determined by the 
sequential least squares algorithm. The multivari- 


able algorithm can be obtained by replacing scalar 
quantities with their corresponding matrices and 
by following the steps outlined in (2) and (3) 
(Ljung [8]). In particular, the estimate and the 
data matrices take the following forms: 

fl-Mi A 2 ••• A na B 0 B, ••• (12) 

n = [->■(*- n 2) ••• 

-y(t - na) u(t - 1) u(t - 2) 

X ■ • • u(t - nb)] r . ( 13) 

2.3. Identification of interface shape by image proc- 
essing 

When the material (PbBr : ) solidifies, there is a 
distinct color change from dark red (liquid) to a 
light yellow (solid). The image processing meth- 
odology makes use of this information to deter- 
mine the shape of the interface. For acousto-optic 
materials, such as lead halides, there is also an 
observable change in color that can be used to 
detect the interface. 

2.3.1. Detection of interface 

Hue, saturation and intensity 7 (HSI) are the 
three basic components that define the color. 
Significant changes in any one of these compo- 
nents can be used to detect the edges of an 
object. There exist many different methodologies 
to accomplish this. These may include (Gonzales 
and Wintz [9] and Wahl [10]): 

(1) histogram based segmentation; 

(2) edge detection by gradient; 

(3) edge models. 

In this study, we investigated all three methods 
of edge detection. We selected the edge model 
based detection procedure for the interface quan- 
tification because of its robustness. 

2.3.2. Heuckel edge model 

Edge models are a means of fitting parameter- 
ized templates to an image. They offer an ap- 
proach to position a particular optimal cun ; e to 
fit a certain bounded region within an image. The 
method used was originally known as the Heuckel 
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Fig. 1. Flow chan for image processing. 


edge model [11]. The model is used to find an 
optimal line within a region and it gives an indi- 
cation of the strength of the edge. The method 
we used is a two-step process as shown in fig. 1. 
The steps are described below. 


Step I: Optimum edge angle 

This method was derived using Heuckel's edge 
model by Merro and Vassey [12]. They show that 
the edge angle of the best fit line (a) of an 8 X 8 
region can be found by convolving two orthonor- 
mal bases, h ^ and h 2 , as shown below': 


" EE/a, *) +) 


(14) 


coordinates within the operator's window', and a 
is the edge angle with respect to h v The orthogo- 
nal bases, h } and h 2 , are shown in fig. 1. 

Step 2: Position of the edge 

Now that the optimal edge angle (a) is known, 
it is necessary to find the position of the edge. 
The angle a is used to make an edge model 
template g(f, tp). Next, the squared error be- 
tween the image f (£ . <//) and the template g{£, ip) 
is minimized by the following error measure: 

* = EE [/a, ~sa, ^)]v 

+g(e, </')*] • 0 6 ) 

This measure can be minimized by maximizing 

p * EE/a> *i>)- O") 

The point (£, ip) that maximizes (17) determines 
the edge position. Details of the algorithm can be 
found in Sharpless et al. [13]. 

2.3.3. Determination of interface shape 

Once the solid-liquid interface is found, it 
becomes necessary to transform the interface in- 
formation into a more directly useful form. This 
transform should compress the data into one or 
more control variables which may be used to 
determine the quality of the crystal being pro- 
duced. This shape may fall into one of the follow'- 
ing three classes: flat, concave, or convex. 

The shape can be determined by using radius 
of curvature at several points along the edge. The 
sign of the radius of curvature gives a measure- 
ment of convexity or concavity. The magnitude of 
the radius is a gauge of the amount of curvature. 

3. Furnace control 

3.7. Seif-tuned pole placement temperature con- 
troller ; 5750 case 

The control law for the process is assumed to 
be 

U(r)-T(s~ l ) r(,)-S(z-') y(t ). 

(18) 

where T(z ~ '), and S(z _1 ) are polvnomi- 


where f((, ip) is the image intensity' of the local 
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als to be determined. Substituting (18) into the 
process dynamics (1) and solving for y(r), gives 

A(z-')R{z-')+B(z-')S(z-') nt) - 

(19) 

If the desired characteristic polynomial for the 
feedback control system is P(z~ l \ it follows from 
(19) that 

A{z-') R(z-')+B(z-') S(z - 1 ) 

-F(z~') P(z - ] ), (20) 

where 7(z _1 ) is an observer polynomial. To 
guarantee zero steady state error for step input, 
the controller is designed to have an integrator, 
i.e., 

R,{z-'). (21) 

Then, by defining proper F(z~ T ) and P(z _1 ), 
R { (z~ : ) and 5(z~ ! ) can be solved by equating 
corresponding powers of z -1 in (20). 

The reference filter 7(z~ 1 ) is solved by con- 
sidering the steady state relation between r(f) 
and y(r). Since R( 1) = 0, it follows from (19) that 
in the steady state: 

, 7X1) ■ 

' V B(l)S(\) r S( l) r ‘ ( ' } 

If 7(1) = 5(1), then y = r in the steady state. 

3.2. Self-tuned pole placement temperature con- 
troller, MIMO case 

The control law used in the multivariable case 
has the same form as in (18). However, the scalar 
polynomials R, 7, and 5 are replaced by matrix 
polynomials. Substituting (18) into (6). the closed 
loop process equation becomes 

■Ay(t) = BR~ } [Tr(t) — 5y(r)] +n(r). (23) 

Arranging for y(r) gives 

[A +BR-' t S} y(r) = z _1 BR~'T r(i) + n{t), 

(24) 


R-'S = SR~\ (25) 

with nr = nr and ns = ns. Then, (24) can be writ- 

ten as 

[ AR +55]^" ] y(f) = z~ l BR~ ] T r(t) + n(t). 

(26) 

The solution for R, 5, and T follows a similar 
approach as in the case of S1SO design, by setting 

[.4f? + BS] =/ + /V-' +P,z- ] + -r P np z-' :p 

= / 5 ( 2 ' 1 ), ( 2 ") 

where P(z~ ] ) is the designer's closed loop char- 
acteristic polynomial with each P 2 as an 8 X 8 
matrix. The order of P(z~ l ) is determined by the 
model order and the controller order as 

np = max{ na 4- nr , 1 + nb + ns) . (28) 

By imposing zero steady state error, the 7 
polynomial can be solved similar to the SISO 
case. By setting 

R = (I-Iz-')R { , (29) 

the steady state solution for (26) is 

[,4(1) R(l)+B(l) 5(1)] R-\ 1) y(r) 

-fi(l) fr’(l) T(l) f(t). (30) 

Using (25) and (29), (30) becomes 

tf' ] (l) 5(1) y = R~\ 1) 7(1) r . (31) 

With the choice of 7(1) = 5(1), eq. (31) becomes 
y(t) = r(t) in the steady state. 

In order to implement the control law, R and 
5 must be transformed to R and 5 using (25), i.e. 

SR = RS. (32) 

4. Experimental results 

The major motivations behind the experimen- 
tal study are: 

(A) Comparison between the SISO and MIMO 
self tuning temperature control algorithms. 




Fig. 2. Temperature control hardware: (1) computer controller with A/D and D/A boards; (2) thermocouple transducers; (3) duty 
cycle based heater drivers; (4) translational mechanism controller; (5) independent data acquisition system; (6) independent on-off 
controller for safety’; (7) transparent furnace; (8) thermocouples; (9) heater connections; (10) translational mechanism. 


(B) Image processing based quantification of the and in each zone the temperature is measured by 
interface. three K type thermocouples. Based on 12 bit 

Fig. 2 shows the temperature control system, A/D conversion and a measurement range of 
The furnace has eight interacting heating zones 700 C C the temperature reading resolution is 



Fig. 3. Image processing hardware: (1) computer; (2) intermittent video recorder and video monitor; (3) TV 7 camera and front-end 

lens system: (4) furnace. 
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Fig. 4. Furnace operating under S1SO control. 


MIMO CONTROL 



TIME (sec) 

Fig. 5. Furnace operating under MIMO control. 
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Fig. 6. Color component contour plots. 



0.12°C. The heater drivers are zero crossing, duty 
cycle modulated solid state switches. Each can be 
driven from a D/A port of the computer with a 
12 bit resolution. The main disturbance on the 
control system is the motion of the ampoule. This 
is provided by a computer controlled translational 
mechanism. The motion resolution on this system 
is approximately 1 per step. 

The imaging system is shown in fig. 3. It can 
provide a view of the interface with an eight-bit 
resolution on each color. The total number of 
pixels representing the image is 512 X 480. 


Figs. 4 and 5 show the performance of the 
control system both in SISO and in MIMO self- 
tuning mode. In each case, two different experi- 
ments are performed: 

(1) For the first experiment, all reference tem- 
peratures are first decreased and then increased 
simultaneously by 50°C. The aim of this experi- 
ment is to see if the response follows the curve 
imposed by the desired pole locations. 

(2) The second experiment is designed to see the 
effect of zone to zone interactions. Only one 
reference temperature is changed by 75°C, all 



Fig. 7. End result of interface identification. 
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others are kept constant. For both SISO and 
MIMO control, the desired poles are placed at 
location z = 0.9. This choice results in an over- 
damped response. 

The results of these experiments lead to the 
following conclusions: 

(a) In MIMO control, all zones behave in a simi- 
lar pattern. This pattern is imposed by the de- 
sired pole locations. In the case of SISO control, 
zones 6, 7, and 8 show underdamped responses. 
This illustrates that the SISO pole placement 
design cannot handle zone-to-zone interaction as 
expected. 

(b) In the case of MIMO control, if one zone's 
reference temperature is changed, this change 
has an effect on the neighboring zones. However, 
the effect is not long lasting in comparison with 
the same changes introduced in the case of SISO 
control. 

(c) There is no noticeable difference between the 
SISO and MIMO controllers in the steady state 
operation. 

The image quantification starts with the red, 
green, and blue components of the image. Fig. 6 
shows these components around the interface. As 
it can be observed, the green component provides 
the greatest contrast around the interface. This is 
illustrated by the close fit contour lines which 
show the interface better in the green than in the 
red and blue components. The Heuckel algorithm 
outlined in fig. 1 is applied to the green compo- 
nent of the picture. The resulting image is shown 
in fig. 7 where the identified interface is marked 
by a bold black line. The speed of the algorithm 
is approximately 5 s, which is properly suited to 
the time step within the temperature control al- 
gorithm. During this time the algorithm finds the 
interface within an 8 x 255 pixel region. 

5. Conclusion 

This paper presents an innovative approach, 
employing image processing techniques in con- 
nection with the conventional identification pro- 
cedures to identify and control the dynamics of a 
crystal while it is growing inside a transparent 
furnace. 


A multizone transparent crystal growth fur- 
nace is controlled by an adaptive pole placement 
controller. Changes in the furnace dynamics and 
the interaction between neighboring heating 
zones are estimated by an on-line least squares 
algorithm. The pole placement control algorithm 
is extensively tested on the growth of lead bro- 
mide crystals. 

For an explicit model of the heat interactions 
between heating zones, a multi-input-multi-out- 
put (MIMO) model is determined. A self-tuning 
multivariable regulator is developed and applied 
to the furnace. Comparison between the resulting 
single-input-single-output (SISO) and MIMO 
controllers is performed. This comparison shows 
that the MIMO controller provides better control 
in eliminating zone-to-zone heat interactions. This 
is important because it avoids local overheating 
of any zone due to interaction. Thereby, a much 
smoother temperature profile structure can be 
imposed. 

Several different image processing routines 
have been applied to find the solid-liquid inter- 
face while a crystal is growing within a transpar- 
ent furnace. A technique based on the Heuckel 
edge model is found to be the most effective 
means to determine the position and shape of the 
interface. This shape information can be used to 
assess the qualify of a crystal while it is growing. 
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In our letter dated, May 26, 1992, we requested certain documents in order to 
closeout the subject cooperative agreement and retire our file. According to 
our records we have not received the following materials: 


1. Delivery of the Final Technical Report distributed as follows: 

a. Orginial to Grants Officer, Mail Stop 500-215. 

b. Copy to NASA Technical Officer. 

c. Copy to NASA Center for Aerospace Information, P.0. Box 8757, BWI 
Airport, MD 21240. 

2. Submission of the final equipment inventory report, a negative report is 
required should no inventory exist. 

If complaince has been made, please send copies to the attention of Pat 
Justice, Procurement Assistant, Mail Stop 500-315. Thank you for your 
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